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Abstract. Moving puncture simulations of black hole binaries rely on a specific gauge 
choice that leads to approximately stationary coordinates near each black hole. Part 
of the shift condition is a damping parameter, which has to be properly chosen for 
stable evolutions. However, a constant damping parameter does not account for the 
difference in mass in unequal mass binaries. We introduce a position dependent shift 
damping that addresses this problem. Although the coordinates change, the changes 
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in the extracted gravitational waves are small 

U . 
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*0 . 1. Introduction 

Early numerical relativity simulations using a 3+1 split of the Einstein equations suffered 
from so-called slice stretching, an effect which occurs when using singularity avoiding 
; slicing together with a vanishing shift. The slices become highly distorted when time 
marches on in the outer regions of the grid but slows down in the vicinity of the black 
hole. It became clear that a non-vanishing, outwards pointing shift vector would be 
required in order to redistribute grid points and also to prevent grid points from falling 
into the black hole. Inspired by Balakrishna et al. pQ, Alcubierre et al. [21 [3] combined 
the 1+log slicing condition with a dynamical shift condition called gamma-driver. These 
gauge conditions successfully prevented slice stretching in black hole simulations using 
excision. It turned out that such gauge conditions could be used also for fixed punctures 
with slight modifications to keep the puncture from evolving [I]. The fixed-puncture 
modification was removed in [51 E] when the moving puncture method was introduced. 
1+log slicing with gamma-driver shift succeeds in moving the puncture freely through 
the grid while simultaneously avoiding slice-stretching. The basic reason for the success 
of this gauge condition is that when the slices start to stretch, the shift vector counteracts 
by pulling out grid points from the region near the black hole. 

In this paper we focus on the dissipation or damping parameter in the gamma- 
driver shift condition, which plays an important role in the success of this gauge. In 
order to reduce oscillations in the shift vector, the authors of [3] noticed the necessity 
of a damping term in the shift condition. Adjusting the strength of the damping via a 
damping parameter was found to allow freezing of the evolution at late times in [4] and to 
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avoid drifts in metric variables in j7] . Additionally, the value of the damping coefficient 
was found to affect the coordinate location of the apparent horizon and therefore the 
resolution of the black hole on the numerical grid [3 [8] . The right choice of the damping 
value is therefore important if one wants to resolve the black hole properly while still 
driving the coordinates to a frame where they are stationary when the physical situation 
is stationary and hence obtain a stable evolution. 

The specific value of the damping parameter has to be adapted to the black hole 
mass in order to obtain long term stable evolutions. If the damping parameter is 
either too small or too large, there are unwanted oscillations or a coordinate instability, 
respectively. In binary simulations, a typical choice is a constant value of roughly 2/M, 
where M is the total mass of the system. However, using a constant damping parameter 
for black hole binaries with unequal masses leads to a fundamental problem. With a 
constant damping parameter, the effective damping near each black hole is asymmetric 
for unequal black hole masses since the damping parameter has dimensions 1/M. For 
large mass ratios, this asymmetry in the grid can be so large that simulations fail because 
the damping may become too large for one of the black holes. This is one of the reasons 
why the highest mass ratio that has been successfully simulated up to now is 10 : 1 [9]. 

Advantageous would be a position-dependent damping parameter that adapts to 
the local mass, in particular such that in the vicinity of the i^ 1 puncture with mass 
Mi its value approaches 1/Mj. It was noticed before [HE] that a damping coefficient 
adapted to the various parameters of the simulation would be beneficial. In [TO] a 
position-dependent formula was introduced for head-on collisions of black holes, which 
to our knowledge was only used in one other publication [11], prior to the moving 
puncture framework. In this paper, we take first steps towards a position-dependent 
damping parameter for moving punctures. As a consequence, the local coordinates 
change compared to standard simulations, but this does not significantly affect gauge 
invariant quantities like the extracted waves as we discuss below. 

2. Dynamical damping in the shift equation 

2.1. Numerical setup 

We focus on the gauge condition used in the 3+1 splitting of the Einstein equations, in 
particular on the condition for the shift vector. 

The slices are determined by the 1+log slicing condition [12] for the lapse function 

a. 



where K is the trace of the extrinsic curvature. The coordinates of a given slice are 
governed by the gamma-driver shift condition introduced in pE] as 



d a = —2aK, 



(1) 




h r - Tjsdo?, 



(2) 
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where P are the contracted Christoffel symbols of the conformal metric 7^, j3 l is the 
shift vector and r] s is the damping coefficient we will discuss in this publication. In 
Eqs. (JTJ and (TSJ), d is defined as d = d t — (3 % di as suggested by [UH HH EE]- 

Examining the physical dimensions, we see that [f3 l ] = 1 and [do] = 1/M, where M 
is the mass (e.g. the total mass of the spacetime under consideration). For this reason, 
the second term on the right hand side of equation (T5]) requires the damping parameter 
to carry units, 

M = ±. (3) 

In simulations of a single Schwarzschild puncture of mass Mi, we typically choose 
a damping parameter of r] s ~ 1 /Mi for obtaining enough damping in the shift without 
producing instabilities. In numerical experiments for a Schwarzschild puncture (to be 
discussed elsewhere), we find that < r] s ;$ 3.5 /M 1 is necessary for a stable and 
convergent numerical evolution. Some minimal amount of damping is important to 
suppress noise in the gauge when a puncture is moving. On the other hand, if i] s is too 
large then there are gauge instabilities, leading to a loss of convergence and to instability 
of the entire numerical evolution. Furthermore, early simulations for fixed punctures 
also found that n s should take values around 1/M, where M is the total mass, to avoid 
long-term coordinate drifts at the outer boundary [I]. 

In simulations of black hole binaries with total mass M = Mi + M 2 , we usually set 
r] s = 2/M which has been found to work well in equal mass binaries simulations. (For 
equal masses, Mi = M 2 = M/2, so near one of the punctures the value of i] s discussed 
above for Schwarzschild becomes r] s = I /Mi = 2/M.) For unequal mass binaries, the 
different black holes tolerate different ranges of r] s according to the above statement 
about single black holes. Ideally, r] s should be ~ 1/Mj, which cannot be accomplished 
simultaneously for unequal masses using a constant value of rj s = 2/M. In fact, for the 
mass ratio 1:10 in [9], the choice r] s = 2/M failed, but a smaller value for n s was chosen 
such that i] s ;$ 3.5/ Mi for both i = 1 and i = 2. 

To overcome the conflicts between punctures with different masses in evolutions of 
two or more black holes, we suggest to construct a non-constant, position-dependent 
damping parameter which knows about the position and mass of each puncture and 
takes a suitable value at every grid point. 

2.2. Using ip~ 2 to determine the position of the punctures 

We thus desire a definition of T] s which respects the unit requirements found in Eq. (J3]) 
and which asymptotes to specifiable values at the location of the punctures and at 
infinity. Typical values are r/ s = 1/M; at the z^ 1 black hole and n s = 2/M at large 
distances. This can be achieved by determining n s through a position dependent function 
defined on the whole grid instead of using a constant as before. We desire a smooth r] s 
which avoids modes which travel at superluminal speeds. Since we use the Baumgarte- 
Shapiro-Shibata-Nakamura (BSSN) system of Einstein's equations [TBI D2]> we want 
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the form of 7] s to depend only on the BSSN variables in a way that does not change the 
principal part of the differential operators. 

In this paper, we choose to use the conformal factor if), which contains information 
about the locations and masses of the punctures. The formula we will use for determining 
the damping coefficient r] s (f) is 

with 7* J the inverse of the conformal 3-metric and Rq a dimensionless constant. While 
ip, an d Rq are dimensionless, the partial derivative introduces the appropriate 
dependence on the mass since [di\ = 1/M and hence [7/ 8 (f)] = 1/M. 

For a single Schwarzschild puncture of mass M located at r = the behavior of 
Eq. (jl]) near the puncture and near infinity is as follows. According to [18], for small 
radii r (near the puncture) the conformal factor asymptotically equals 

ip~ 2 ~ pir (5) 



w) = Ro^-i , (4) 



for a known constant p\. The next to leading order behavior is less simple [19]. The 
point r = corresponds to a sphere with finite areal radius Ro, 

R = lim^V = — = R M. (6) 

Numerically, Rq ~ 1.31. The inverse of the conformal metric behaves like 

f j ~ 5 ij . (7) 
Therefore, we find for small r 

^/Y j dii/j~ 2 djip~ 2 ~ pi = -J 1 — (8) 

RqM 

and 

(1 -?T 2 ) 2 ~ (1 -pir) 2 ~ 1 (9) 

when keeping only leading order terms in r. Equations (jSj) and combine according 
to (jl]) to give 

r) s ( r = 0) = 1/M. (10) 
For large r we can expand the conformal factor in powers of 1/r, 



resulting in 



and 



2r J 



y/fWi^djl/j- 2 = (12) 



b M/r 2 Ro 

Vs(r - oo) - R 0WM - 2 = -. (13) 
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In summary, Eq. (J4J) leads to 



for a single puncture at r = 0. Note that using Eq. (@J in Eq. (TSJ) does not affect the 
principal part of (j2J). Therefore, the system remains strongly hyperbolic, same as for 
T] s = const, according to [El [13] . 

3. Results 

Our Eq. (j4]) analytically gives the desired 1/M behavior near the puncture and near 
infinity for a single, non-spinning and non-moving puncture. Now it remains to be 
tested whether these properties persist in actual numerical simulations, especially for 
unequal mass binaries. 

Simulations are performed with the BAM code described in [7J [20]. The code 
uses the BSSN formulation of Einstein's equations and employs the moving puncture 
framework [5J E]. Spatial derivatives are 6*^ order accurate and time integration is 
performed using the 4^ order Runge-Kutta scheme. The numerical grid is composed 
of nested boxes with increasing resolution, where the boxes of highest resolution are 
centered around the black holes. These boxes are advanced in time with Berger-Oligar 
time stepping |21j. We are using puncture initial data with Bowen-York extrinsic 
curvature and solve the Hamiltonian constraint using a pseudospectral collocation 
method described in [22]. The momentum parameter in the Bowen-York extrinsic 
curvature is chosen such that we obtain quasi-circular orbits in our binary simulations 
using the method of [23] . 

For binary simulations with unequal masses, we will use the mass ratio q = M2/M1 
to denote the runs, Mj being the bare mass of the i^ 1 puncture. The physical masses 
of the punctures (obtained after solving the constraints) differ by less than 10% from 
the bare masses for the orbits considered here, so the t] s values derived for a single 
puncture should remain valid. When comparing simulations run with rj s = 2/M and 
7] s (r) following Eq. (jlj) we will refer to them as "standard" and "new" or "dynamical" 
gauge, respectively, throughout this paper. 

3.1. Single Schwarzschild Black Hole 

In order to test the 1/M-behavior of (jl]) near the puncture and infinity, we first 
performed a series of evolutions for a time of 100 M of a single, non-spinning puncture 
while varying its mass. We then measured the value of r] s (r) near the puncture and 
at the outer boundary of the grid and compared these values to the limits (TBI) . The 
data points in Fig. [1] correspond to these measurements while the lines are fits to the 
numerical data. The values of i] s (r) near the puncture as a function of total mass M 
are fitted to f] s {M) = 1.05/M which agrees well with the analytical limit r — > of (THj) . 
Fitting to rjs measured near the physical boundary of the grid reveals rj s (M) = 1.311/M 




(14) 
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Figure 1. Numerical test of the analytical limits IT4|) of rj s (r) using single, non- 
spinning punctures with different masses. Shown are the values of ri s (r) near the 
puncture (gray squares) and at the outer boundary (black dots). Note that here M is 
identified with a dimcnsionlcss number, so w s {r) is dimensionless as well. The fits to the 
data points are consistent with the analytic prediction. Numerically, t] s (M) — 1.05/M 
(gray line) near the puncture and r) 8 {M) = 1.311/M (black line) at the outer boundary. 
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Figure 2. ^-component of the shift vector (upper panel) and xx-componcnt of the 
conformal 3-metric (lower panel) of a single, non-spinning puncture at time t = 100 M, 
where the simulations have reached a stationary state. The dashed black curves use 
dynamical damping, Eq. (Tj|, the gray ones use rj s = 2.0/M in the shift condition 
Eq. ©. 



and therefore fulfills the limit r — > oo of IHM even though the outer boundary is situated 
only at 130 M. 

Using a modified shift condition, the shift itself will, of course, change. We compare 
the x-component of the shift vector for using t] s = 2.0/M and n s (r) in Fig. [2]. A change 
in the shift implies a change of the coordinates and therefore, coordinate dependent 
quantities will change, too. As an example, the xx-component of the conformal 3-metric, 
j xx , is compared for rj s = 2/M and r] s (f) in the lower panel of Fig. [2l The comparisons 
are made at time t = 100 M, when the simulations have reached a stationary state. 

The changes in the shift should only affect the coordinates and coordinate 
independent quantities should not change. This can be examined by looking at a scalar 
as a function of another scalar, e.g. the lapse a as a function of extrinsic curvature K, 
a = a(K). Both scalars should see the same coordinate drifts and therefore, no changes 
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Figure 3. The lapse function a as a function of extrinsic curvature K for a single, 
non-spinning and non-moving puncture after a time t = 50 .M. We compare using r; s (r) 
(black, dashed line) and y s = 2.0/M (gay line). The two curves lie perfectly on top 
of each other and are therefore indistinguishable. The insets show lapse (upper panel) 
and extrinsic curvature (lower panel) as functions of distance from the puncture. 

are expected in a(K). Figure [3] confirms this expectation. The two curves a(K) for 
7] s = 2.0/M and ri s (r) are lying perfectly on top of each other. We therefore believe 
that using the dynamical damping introduces only coordinate changes in our puncture 
simulations. 

3.2. Black hole binary with equal masses 

While Eq. ((!]) has been introduced in order to allow for numerical simulations of two 
black holes with highly different masses, we first apply it to equal mass simulations in 
order to perform several consistency checks. 

The (first order) coordinate independent quantity to look at in binary simulations 
is the Newman-Penrose scalar \l/ 4 . We use ^ 4 for the extraction of gravitational waves 
(see [7] for details of the wave extraction algorithm), decomposed into modes using 
spin- weighted spherical harmonics ■ Since is only first-order gauge invariant and 
we furthermore extract waves at a finite, fixed coordinate radius, it is a priori an open 
question how much the changes in the shift affect the wave forms. 

As the most dominant mode of ^4 in an equal mass simulation is the I = \m\ = 2 
mode, its real part multiplied by the extraction radius (rex = 90 M in this case) is 
displayed in Fig. HI We look at amplitude and phase of this mode in Figs. |5] and [6j The 
initial separation was chosen to be D = 7 M. The black holes complete about 3 orbits. 
Three different resolutions are used corresponding to the three different colors in Figs. HI 
|5]and[6l We use the number of grid points in the inner boxes (centered around the black 
holes) to denote the different resolutions. The grid configurations, in the terminology 
of [7J, are 0[5 x 56 : 5 x 112 : 6], 0[5 x 64 : 5 x 128 : 6], and 0[5 x 72 : 5 x 144 : 6] 
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Figure 4. Real part of the 22-mode of #4 times extraction radius rex for an equal 
mass binary with initial separation D = 7 M using r/ s = 2.0/M (solid lines) and rj s (r) 
following Eq. (Tj[]) (dashed lines) in three different resolutions (blue, red, green lines) 
according to the grid configurations described in the text. 

which corresponds to resolutions on the finest grids of 3M/112, 3M/128 and M/48, 
respectively. These are the grid configurations used in [20] . 

In Fig. [5], we compare the amplitude A 2 2 in the standard gauge, rj s = 2.0/M, 
displayed as solid lines, to the new one, Eq. PJ, plotted as dashed lines. We find that 
the differences between standard and new gauge for a given grid resolution are much 
smaller than differences due to using different resolutions. This strengthens the belief 
that we only introduced coordinate changes to the system when using Eq. (jlj). The 
maximum relative deviation between the amplitudes A22 of old and new gauge amounts 
to about 3% for the lowest resolution (N = 56) and decreases with increasing resolution, 
which can be seen in the inset of Fig. [5j For the phase, the absolute differences are not 
visible by eye and therefore, we only plot the relative deviations between $22 in the 
standard and new gauge for the three resolutions. For the lowest resolution (N = 56), 
the maximum deviation is only 0.35%. As for the amplitude, this deviation decreases 
with increases in resolution. This shows how the differences in the waveforms disappear 
with increasing resolution. 

The fact that there actually are differences visible in the waves, though very small 
ones, is not surprising when considering the way we extract gravitational waves. We fix 
a certain extraction radius and compute the Newman-Penrose scalar on a sphere of this 
radius. The radius itself is coordinate dependent and we are comparing \l/ 4 extracted 
at slightly different radii in the standard and new gauges. In future work we plan to 
compare wave forms extrapolated in radius to infinity, although it is worth noting how 
small the deviations are without additional processing. 

3.3. Black hole binary with mass ratio J^:l 

After having examined the influence of using a dynamical damping coefficient r\ s {f) for 
an equal mass binary, the next step is to look at its behavior for unequal masses. The 
following results are obtained from a simulation of two black holes with mass ratio q = 4 
and initial separation D = 7 M. We used the grid configurations <p[5 x N : 7 x 2N : 6] 
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Figure 5. Amplitude of the 22-mode of '5 4 for the same binary as in Fig. ([4]) using 
rj s = 2.0/M (solid lines) and rj s (f) (dashed lines) in three different resolutions (blue, 
red, green lines) according to the grid configurations described in the text. The inset 
shows the relative deviation Aj4 2 2/^22 between the amplitude in the standard and in 
the new gauge, again for the three different resolutions. 
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Figure 6. Relative phase difference of the 22-mode of ^4 for the same binary as in 
Fig. ((4| using r] s — 2/M and y s (r)- Compared are three different resolutions (blue, red 
and green lines) according to the grid configurations described in the text. 

with N = 72,80 which have also been used in [24] for mass ratio 4:1. An interesting 
question in this context is how Eq. (J4J) behaves for a simulation with two punctures 
with different masses. The analytical behavior (fT4"|) was deduced for a single, non- 
moving, stationary puncture but now we are using it for two moving punctures, which 
are during most of the simulations far from having reached a stationary state globally, 
but approximately stationary locally at the punctures. 

Figure [7J illustrates the distribution of r] s (r) between the two punctures. For 
convenience, the conformal factor = Inip is also plotted in order to indicate the 
positions of the punctures via its maxima (the divergences are not resolved). The 
snapshot is taken at a time during the simulation when the punctures are still well 
separated. Similar to the simulations of a single puncture, according to ( 1T4")) we expect 
to find r) s (r) ~ 2 near the puncture with mass Mi = 0.5 and i] s (r) ~ 0.5 in the vicinity 
of the second puncture with M 2 = 2.0. Near the outer boundary, i] s is supposed to 
take the value 1.31/(Mi + M2) = 0.52. Here the Mj are chosen to be dimensionless, so 
r] s (r) is dimensionless as well. Figure [7] confirms that we do obtain the expected values, 
although they are not reached exactly. The latter is not a problem as simulations work 
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Figure 7. rj s (r) (black line) for two punctures of bare masses Mi = 0.5 and Mi = 2.0 
and the conformal factor <fi (gray line) whose maxima show the current positions of 
the punctures. In this plot the Mi are dimensionless, so rj s (r) is dimensionless as 
well. The snapshot is taken right after the beginning of the q — 4 simulation with 
initial separation D — 5 M described in the text. The smaller black hole is at position 
d = —3.6 M, the larger one at position d — 1.2 M. 
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Figure 8. Real part of the 22-mode of '5 4 multiplied by the extraction radius rex 
for g = 4 and initial separation D = 5M runs. Compared are results for employing 
j] s = 2.0/M (solid lines) and rj s (r) (dashed lines) in two different resolutions (red and 
green lines) according to the grid configurations described in the text. 

nicely as long as r] s (f) is in the right range for each black hole. For this reason, Eq. (j4j) 
also seems to work rather nicely for two punctures with unequal masses. 

As we did in the equal mass case in section 13.21 we compare the 22-mode of ^4 
in the new gauge with the standard gauge. To see how small the differences using rj s 
or i] s (r) are, we plot its real part using two different resolutions which correspond to 
the two different colors in Fig. [HJ Figure [9] shows the amplitude for the two different 
resolutions. The inset gives the relative differences between amplitudes in the standard 
and new gauge. The maximum relative deviation appears for the lower resolution and 
amounts to about 3%. The high resolution gives 0.5% relative difference. The phases 
in standard and new gauge are compared in Fig. [101 Again, we show only the relative 
deviations as the absolute ones are too small to be seen. We find relative differences of 
up to 0.4% for the lower resolution and only 0.1% for the high one. This confirms that 
we are changing only the coordinates, as we found before in Sections 13. II and 13.21 

While the invariance of the waveforms is the most important feature of the new 
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Figure 9. Amplitude of the 22-mode of ^4 for the same runs as in Fig. [5] We 
compare results using r/ s — 2.0/M (solid lines) and rj s (f) (dashed lines) in two different 
resolutions (red and green lines) according to the grid configurations described in the 
text. The inset shows the relative deviation AA22/^22 between the amplitude in the 
standard and in the new gauge, again for the same two resolutions. 
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Figure 10. Relative phase difference of the 22-mode of ^4 for the same run as in 
Fig-El?7s = 2.0/M and rj s (r). Compared are two different resolutions according to the 
grid configurations described in the text. 

gauge rj s (r), it is illuminating to examine how the black holes are represented on the 
numerical grid. To this end, the apparent horizons (AH) are computed for both gauges 
in Fig. [TTJ We show the result in the (x, y)-plane, in which the orbital plane lies. For 
clarity, the slices through the apparent horizons are only shown at 4 different times. In 
the beginning of the simulations, the AH pertaining to the same black hole are lying 
on top of each other. With time, they separate as the coordinates become more and 
more different in the two simulations. Two observations can be made. First, the ratio 
between the coordinate area of the AH of the larger black hole and the one of the smaller 
black hole is larger in the simulation using rj s = 2.0/M. This means the black holes are 
represented more equally on the grid in the simulation using Eq. (j3j). This fact can be 
seen even more clearly in Fig. [12] where we plot the coordinate area of the apparent 
horizons comparing the standard gauge (red lines) and the new one (black lines). While 
the coordinate sizes of the smaller black hole (dashed lines) are nearly equal in both 
gauges, the sizes of the larger black hole (solid lines) differ by roughly 2 M 2 . Second, 
the shape of the horizon of the smaller black hole is more and more distorted in the 
rjs = 2.0/M-simulation when the black holes come closer together. This deformation is 
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Figure 11. Comparison of apparent horizons in the orbital plane using i] s — 2.0 /M 
(red lines) and rjs{r) (black lines) at different times during the evolution for a q = 4 
run with initial separation D — 5 M. 




Figure 12. Comparison of the coordinate area of the apparent horizons using 
r/ s = 2.0/M (red lines) and f] s (r) (black lines) over evolution time for a q = 4 run with 
initial separation D = 5 M until shortly before a common apparent horizon appears. 
The dashed lines belong to the smaller black hole whereas the solid lines represent the 
larger black hole. 



not visible in the new coordinates. The progressive stretching of the apparent horizon 
shape and therefore the distortion of the coordinates near the black holes can be a source 
of instabilities, e.g. [25]. Using Eq. (j3J) seems to be profitable in this regard. 

3.4- Behavior ofr] s {f) and influence on the shift vector 

Despite the encouraging results we have seen so far, there is a non-negligible concern 
using Eq. (j4j) in the gamma- driver condition ([2j). Although we do not determine 
the damping coefficient via a wave equation, we see wavy features in rj s (r) traveling 
outwards. These distortions even leave remnants on the grid, especially when they pass 
through a refinement boundary. The form of rj s (x) after different evolution times can 



Toward a dynamical shift condition for unequal mass black hole binary simulations 13 



2.5 
2 

1.5 

l 

0.5 

% 

Figure 13. Form of r] s (f) in re-direction at different times during an equal mass 
binary simulation. Noise travels outwards and leaves strong distortions on the grid. 
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Figure 14. x-componcnt of the shift vector in x-direction (solid curves) and y- 
component of the shift vector in y-dircction (dashed curves) after the merger of two 
equal mass black holes at time t = 500 M when using either the standard gauge 
•q s = 2.0/M (gray lines) or the dynamical one, rj s {f), (black lines). 

be seen in Fig. [13] for the equal mass binary described in Sec. 13.21 The result is similar 
in the q = 4 simulation and even in the Schwarzschild simulation, an outward traveling 
pulse is present, which however does not leave visible distortions on the grid and the 
relative amplitude of which decreases for higher mass. The effort we made before in 
order to achieve the correct value of r) s (f) near the outer boundary seem to be canceled 
out by the disturbed shape we find now. As the peaks travel to a region of the grid 
where we have no punctures, we might take the point of view that the exact value of 
rj s (r) and therefore the distortions are of no importance for our simulations. Indeed the 
oscillations do not translate to oscillations in the shift vector as one might think. In the 
shift, we find no gauge "waves" related to the ones in rj s (r). Nevertheless, there is an 
unusual behavior. After merger, when going away from the punctures the shift does not 
fall off to zero as fast as it does when using rj s = const, but keeps a shoulder (compare 
Fig. (HI) which might lead to an unphysical and unwanted drift of the coordinate system. 
We are planning to investigate these issues in more detail in the future. 
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4. Discussion 

We presented a new approach to determine the coordinates in slices of spacetime for 
binary black hole simulations where we take the distribution of mass over the grid into 
account. We have shown that our approach of determining the damping parameter in 
the gamma-driver condition dynamically via Eq. (T4|) gives stable evolutions and does not 
significantly change the gravitational waves extracted from binary systems of equal or 
unequal masses. Furthermore, the use of Eq. (jlj) in an unequal mass simulation resulted 
in a more regular shape of the apparent horizon of the smaller black hole as the binary 
merges. The coordinate size of the apparent horizons became more uniform with the new 
damping coefficient which is a first step towards representing and resolving black holes 
with different masses equally and hence removing the large asymmetry which usually 
distorts the numerical grid in unequal mass simulations. We found gauge waves in our 
damping coefficient which might affect the stability in very long-term simulations and 
lead to coordinate drifts after the merger of the binary. We will address these issues in 
a future publication [26] . 
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